scpd14.oc <- octav(scpdSM[14,])
par(mfrow=c(4,4), mar=c(3,3,2,0)+0.1)
for (i in 1:nrow(scpdSM)){
plot(rad(scpdSM[i,]), prop=F, cex.axis=0.5, xlab="", ylab="", xlim=c(1, 13), type="b", pch=16)
mtext(side=2, line=2, "Number of Indiv", cex=0.5)
mtext(side=3, line=0.5, paste0("Time=", scpdAges$age[i]), cex=0.75)
mtext(side=1, line=2, "Rank", cex=0.5)
}
# examine the two endpoints
scpd1.oc <- octav(scpdSM[1,])
scpd14.oc <- octav(scpdSM[14,])
par(mfrow=c(4,4), mar=c(3,3,2,0)+0.1)
for (i in 1:nrow(scpdSM)){
plot(rad(scpdSM[i,]), prop=F, cex=0.5, xlab="", ylab="", xlim=c(1, 13), type="b", pch=16)
mtext(side=2, line=2, "Number of Indiv", cex=0.5)
mtext(side=3, line=0.5, paste0("Time=", scpdAges$age[i]), cex=0.75)
mtext(side=1, line=2, "Rank", cex=0.5)
}
?rad
?fitrad
scpd1.gs <- fitrad(scpdSM[1,], 'gs')
scpdSM[1,]
scpdSM[1,]>0
scpdSM[scpdSM[1,]>0]
scpdSM[1,scpdSM[1,]>0]
scpd1.gs <- fitrad(scpdSM[1,scpdSM[1,]>0], 'gs')
?fitsad
scpd1.bs <- fitsad(scpdSM[1,scpdSM[1,]>0], 'bs')
scpd1.ls <- fitsad(scpdSM[1,scpdSM[1,]>0], 'ls')
scpd1.lnorm <- fitsad(scpdSM[1,scpdSM[1,]>0], 'lnorm')
scpd1.bs.oc <- octavpred(scpd1.bs)
scpd1.bs.oc <- octavpred(scpd1.bs)
scpd1.ls.oc <- octavpred(scpd1.ls)
scpd1.ls.oc <- octavpred(scpd1.ls)
scpd1.lnorm.oc <- octavpred(scpd1.lnorm)
scpd1.oc
plot(scpd1.oc)
lines(scpd1.bs.oc, col="red")
lines(scpd1.ls.oc, col="blue")
lines(scpd1.lnorm.oc, col="green")
legend("topright",
c("Broken stick", "Logseries", "Lognormal"),
lty=1, col=c("red","blue", "green"))
plot(scpd1.oc)
lines(scpd1.bs.oc, col="red")
lines(scpd1.ls.oc, col="blue")
lines(scpd1.lnorm.oc, col="green")
legend("topright",
c("Broken stick", "Logseries", "Lognormal"),
lty=1, col=c("red","blue", "green"))
?AICtab
AICtab(scpd1.bs, scpd1.ls, scpd1.lnorm, base=TRUE)
?fitsad
?octav
# fit the models
scpd1.geom <- fitsad(scpdSM[1,scpdSM[1,]>0], 'geom')
scpd1.ls <- fitsad(scpdSM[1,scpdSM[1,]>0], 'ls')
scpd1.lnorm <- fitsad(scpdSM[1,scpdSM[1,]>0], 'lnorm')
# compare the fit of each model
AICtab(scpd1.geom, scpd1.ls, scpd1.lnorm, base=TRUE)
# transform to appropriate class for plotting
scpd1.geom.oc <- octavpred(scpd1.geom)
scpd1.ls.oc <- octavpred(scpd1.ls)
scpd1.lnorm.oc <- octavpred(scpd1.lnorm)
# plot the three models against one another and the data
plot(scpd1.oc)
lines(scpd1.geom.oc, col="red")
lines(scpd1.ls.oc, col="blue")
lines(scpd1.lnorm.oc, col="green")
legend("topright", c("geometric", "Logseries", "Lognormal"), lty=1, col=c("red","blue", "green"))
ks(ecdf(scpdSM[1,], ecdf(scpdSM[2,])
)
)
?ks.test
ks.test(ecdf(scpdSM[1,]), ecdf(scpdSM[2,]))
ecdf(scpdSM[1,])
str(ecdf(scpdSM[1,]))
summary(ecdf(scpdSM[1,]))
ks.test(ecdf(scpdSM[1,]), ecdf(scpdSM[2,]), alternative="two.sided")
ecdf(scpdSM[1,])
ecdf(scpdSM[1,])$n
ecdf(scpdSM[1,])
summary(ecdf(scpdSM[1,]))
str(ecdf(scpdSM[1,]))
ecdf(scpdSM[1,])$ecdf
ecdf(scpdSM[1,])[1]
ks.test(scpdSM[1,], scpdSM[2,], alternative="two.sided")
?cdf
?ecdf
x<- ecdf(scpdSM[1,]
ks.test(ecdf(scpdSM[1,]), ecdf(scpdSM[2,]), alternative="two.sided")
x
x<- ecdf(scpdSM[1,])
x
x[1:8]
summary(x)
richness
ages
?estimateR
alt_richness
plot(H~ages)
H
H <- diversity(scpdSM, index="shannon")
Simp <- diversity(scpdSM, index="simpson")
Shan_even <- H/log(specnumber(scpdSM))
Simp_even <- Simp/log(specnumber(scpdSM))
plot(H~ages)
plotCols
plotCols <- c("black", "blue",  "green", "red")
plot(H~ages, type="l", ylim=c(0, y.max), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
range(H)
range(Simp)
range(H_even)
range(Shan_even)
range(Simp_even)
(
plot(H~ages, type="l", ylim=c(0, 2), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
lines(Simp~ages, col=plotCols[2])
plotCols <- c("black", "blue",  "green", "red")
plot(H~ages, type="l", ylim=c(0, 2), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
lines(Simp~ages, col=plotCols[2])
lines(Shan_even~ages, col=plotCols[3], lty=2)
lines(Simp_even~ages, col=plotCols[4], lty=2)
plotCols <- c("black", "blue",  "green", "red")
plot(H~ages, type="l", ylim=c(0, 2), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
lines(Simp~ages, col=plotCols[2])
lines(Shan_even~ages, col=plotCols[3], lty=2)
lines(Simp_even~ages, col=plotCols[4], lty=2)
legend("bottomleft", legend=c("Shannon", "Simpson", "Evenness_Shannon", "Evenness_Simpson"), lty=1, col=plotCols, bty="n", cex=0.75)
plotCols <- c("black", "blue",  "green", "red")
plot(H~ages, type="l", ylim=c(0, 2), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
lines(Simp~ages, col=plotCols[2])
lines(Shan_even~ages, col=plotCols[3], lty=2)
lines(Simp_even~ages, col=plotCols[4], lty=2)
legend("bottomleft", legend=c("Shannon", "Simpson", "Evenness_Shannon", "Evenness_Simpson"), lty=c(1,1,2,2), col=plotCols, bty="n", cex=0.75)
pollen.6000
rads(pollen.6000[1,])
rad(pollen.6000[1,])
pollen.6000[1,]
?ra
?rad
rad(pollen.6000[1,])
sad(pollen.6000[1,])
sads(pollen.6000[1,])
octav(pollen.6000[1,])
#uncomment the following line if you haven't installed these packages yet
#install.packages(c('vegan', 'neotoma', 'dplyr', 'sads', 'iNEXT'))
library('vegan')
library('neotoma')
library('dplyr')
library('sads')
library('iNEXT')
samwell_site <- get_site(sitename="Samwell Cave") #there's one site in the database that matches this
samwell_datasets <- get_dataset(samwell_site[[1]]) # this search lists all datasets associated with Samwell Cave. There are 7 different datasets for the site: 4 containing vertebrate faunal information and 3 containing geochronologic information.
# Download info for just the last dataset, 14262.
scpd_dataset <- get_dataset(samwell_site[[1]])$'14262'
scpd_download <- get_download(scpd_dataset)
# Extract the count data and the sample ages
scpdFauna <- counts(scpd_download)[[1]]
scpdAges <- ages(scpd_download)[[1]]
# save as R object for later upload in case no internet or API not working
# save(scpdFauna, scpdAges, file="data/scpdObjects.RData")
# Note: Dataset 14261 contains the underlying radiocarbon dates for the site, in case you want to examine those. You'll need to uncomment the next two lines
# scpd_geochron <- get_geochron(14261)
# scpd_geochron[[1]]$geochron
# examine the final downloaded file
str(scpdFauna) # another possible way to examine the data: glimpse(scpdFauna)
head(scpdFauna)
# What are the rownames showing?
# Let's replace the sample ID in the fauna object (rownames) with the ages
# but only do this if the sample IDs match across objects
if (all(rownames(scpdFauna) == scpdAges$sample.id)){
rownames(scpdFauna) <- ages
}
# set the names of the small mammal genera of interest
# note: this is a clunky workaround - ideally we'd be able to harmonize taxonomy using the database/API tools, but they aren't as good in Neotoma as in PBDB because of the way taxonomy is handled
sm_genus <- c("Scapanus", "Microtus", "Neotoma", "Peromyscus", "Sorex", "Spermophilus", "Tamias", "Tamiasciurus", "Thomomys", "Arborimus", "Reithrodontomys", "Sciurus", "Glaucomys", "Myodes", "Aplodontia")
# create a new dataframe with the merged data
scpdSM <- matrix(nrow=nrow(scpdFauna), ncol=length(sm_genus))
for (i in 1:length(sm_genus)){
cols <- grep(sm_genus[i], colnames(scpdFauna))
if (length(cols)==1){
scpdSM[,i] <- scpdFauna[,cols]
}
if (length(cols) > 1){
scpdSM[,i] <- rowSums(scpdFauna[,cols])
}
}
# clean up dataframe
colnames(scpdSM) <- sm_genus
rownames(scpdSM) <- rownames(scpdFauna)
# we can practice building our own functions to count up the species
countFcn <- function(df){
count <- length(which(df>0))
return(count)
}
raw.S <- apply(scpdSM, 1, countFcn)
# or you can do the same calculation using a built-in function (specnumber) in vegan
specnumber(scpdSM)
range(rowSums(scpdSM))
minN <- min(rowSums(scpdSM))
richness <- rarefy(scpdSM, sample=minN, se=T) # because I specified *se=T*, this function outputs both S and se
S <- richness[1,]
S.se <- richness[2,]
# plot the result through time using base R plotting function
ages <- scpdAges$age
plot(S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, floor(max(S+S.se))+1), xlab="Age (cal BP)", ylab=expression('Richness' %+-% 'se'))
#arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
suppressWarnings(arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90))
# Note: use suppressWarnings with caution.  I added the suppressWarnings tag so that the only thing output is the plot. But I only did this after I confirmed that the warning message was ok. Switch the arrows lines that are commented/uncommented and re-run the plot with that line of code for arrows to view the warning and decide for yourself whether it should be suppressed or not.
plot(raw.S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, max(raw.S)+1), lty=2, col="gray", xlab="Age (cal BP)", ylab='Richness')
lines(S ~ ages, type="l")
arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
# set colors, from cool for the older (LGM) times and red for the warmer (Holocene) times
ageColors <- rainbow(length(ages), start=0, end=4/6)
rarecurve(scpdSM, col=ageColors)
?estimateR # read more about this function, and how ACE and Chao1 are calculated
alt_richness <- estimateR(scpdSM)
# Bind all richness estimates together and compare:
richness_all <- rbind(alt_richness, richness)
# plot them against one another
y.max <- floor(max(richness_all))+1 # set y max for plotting
mainEstimates <- c('S.obs', 'S.chao1', 'S.ACE', 'S')
plotRows <- match(mainEstimates, rownames(richness_all))
plotCols <- c("black", "blue", "green", "red")
i=1
plot(richness_all[plotRows[i],] ~ ages, type="l", ylim=c(0, y.max), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[i])
for (i in 2:length(plotRows)){
lines(richness_all[plotRows[i],] ~ ages, col=plotCols[i])
}
legend("bottomleft", legend=mainEstimates, lty=1, col=plotCols, bty="n", cex=0.75)
# this chunk uses the iNext package
vignette("Introduction", package="iNEXT")
iNextRichness <- iNEXT(t(scpdSM))
# Note 1: we had to transpose the community matrix
# Note 2: this calculates not just richness, but also Shannon and Simpson diversities.
# You can also rarefy within iNEXT using estimateD
D_minN <- estimateD(t(scpdSM), datatype="abundance", base="size", level=NULL)
D_minN
# Note 1: SC = sample coverage, qD is the estimate with lower and upper confidence levels
# Note 2: Richness estimates here are the same as using *rarefy* in vegan. e.g., compare the following:
D_minN[which(D_minN$order == 0),]$qD
richness[1,]
ggiNEXT(iNextRichness, type=1, facet.var="order", color.var="site")
ggiNEXT(iNextRichness, type=1, facet.var="site")
ggiNEXT(iNextRichness, type=2)
vignette("diversity-vegan")
H <- diversity(scpdSM, index="shannon")
Simp <- diversity(scpdSM, index="simpson")
Shan_even <- H/log(specnumber(scpdSM))
Simp_even <- Simp/log(specnumber(scpdSM))
plotCols <- c("black", "blue",  "green", "red")
plot(H~ages, type="l", ylim=c(0, 2), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
lines(Simp~ages, col=plotCols[2])
lines(Shan_even~ages, col=plotCols[3], lty=2)
lines(Simp_even~ages, col=plotCols[4], lty=2)
legend("bottomleft", legend=c("Shannon", "Simpson", "Evenness_Shannon", "Evenness_Simpson"), lty=c(1,1,2,2), col=plotCols, bty="n", cex=0.75)
# load the data and examine the structure
pollen.6000 <- read.csv(file="data/pollen/6000.interp.data_clean_sites.csv")
str(pollen.6000)
# clean up the site/rownames
pollen.6000.sites <- pollen.6000$sites
rownames(pollen.6000) <- pollen.6000$sites
pollen.6000 <- pollen.6000[,-1]
# calculate richness and diversity
S.obs <- specnumber(pollen.6000)
H <- diversity(pollen.6000, "shannon")
Simp <- diversity(pollen.6000, 'simpson')
H_even <- H/log(S.obs)
# get site lat/longs
meta <- read.delim(file="data/pollen/All.site.data-withagemodel-finalv2.txt", sep="\t")
site.longlats <- meta[match(pollen.6000.sites, meta$Handle), c('Longitude', 'Latitude')]
# simple relationships with latitude and longitude
par(mfrow=c(2,2), mar=c(4,4,1,1)+0.1)
plot(H~site.longlats$Latitude, pch=16, xlab="Latitude")
plot(S.obs~site.longlats$Latitude, pch=16, xlab="Latitude")
plot(Simp~site.longlats$Latitude, pch=16, xlab="Latitude")
plot(H_even~site.longlats$Latitude, pch=16, xlab="Latitude")
par(mfrow=c(2,2), mar=c(4,4,1,1)+0.1)
plot(H~site.longlats$Longitude, pch=16, xlab="Longitude")
plot(S.obs~site.longlats$Longitude, pch=16, xlab="Longitude")
plot(Simp~site.longlats$Longitude, pch=16, xlab="Longitude")
plot(H_even~site.longlats$Longitude, pch=16, xlab="Longitude")
# this chunk uses the sads package
vignette('sads_intro')
# examine the two endpoints
scpd1.oc <- octav(scpdSM[1,])
scpd14.oc <- octav(scpdSM[14,])
par(mfrow=c(4,4), mar=c(3,3,2,0)+0.1)
for (i in 1:nrow(scpdSM)){
plot(rad(scpdSM[i,]), prop=F, xlab="", ylab="", xlim=c(1, 13), type="b", pch=16)
mtext(side=2, line=2, "Number of Indiv", cex=0.5)
mtext(side=3, line=0.5, paste0("Time=", ages[i]), cex=0.75)
mtext(side=1, line=2, "Rank", cex=0.5)
}
par(mfrow=c(4,4), mar=c(3,3,2,0)+0.1)
for (i in 1:nrow(scpdSM)){
plot(ecdf(scpdSM[i,]), xlab="", ylab="", xlim=c(1, 13), main="", verticals=T, do.points=F, cex.axis=0.75)
mtext(side=2, line=2, "% of Species", cex=0.5)
mtext(side=3, line=0.5, paste0("Time=", ages[i]), cex=0.75)
mtext(side=1, line=2, "% abundance", cex=0.5)
}
# fit the models
scpd1.geom <- fitsad(scpdSM[1,scpdSM[1,]>0], 'geom')
scpd1.ls <- fitsad(scpdSM[1,scpdSM[1,]>0], 'ls')
scpd1.lnorm <- fitsad(scpdSM[1,scpdSM[1,]>0], 'lnorm')
# compare the fit of each model
AICtab(scpd1.geom, scpd1.ls, scpd1.lnorm, base=TRUE)
# transform to appropriate class for plotting
scpd1.geom.oc <- octavpred(scpd1.geom)
scpd1.ls.oc <- octavpred(scpd1.ls)
scpd1.lnorm.oc <- octavpred(scpd1.lnorm)
# plot the three models against one another and the data
plot(scpd1.oc)
lines(scpd1.geom.oc, col="red")
lines(scpd1.ls.oc, col="blue")
lines(scpd1.lnorm.oc, col="green")
legend("topright", c("geometric", "Logseries", "Lognormal"), lty=1, col=c("red","blue", "green"))
# pollen.datasets <- get_dataset(datasettype="pollen", loc=c(-105, 40, -86, 49), ageyoung = 8000)
# browse(pollen.datasets)
#
# # pollen.downloads <- get_download(pollen.datasets)
# # save(pollen.downloads, file="data/pollen.downloads.RData")
# load("data/pollen.downloads.RData")
#
# # Standardize the taxonomy
# pollen.whit <- compile_taxa(pollen.downloads, list.name="WhitmoreSmall")
# pollen.ages <- ages(pollen.whit)
#
# # find the closest sample to the focal time and save the pollen counts for that sample
# # Note: the loop code not working yet.
# focalTime <- 12000
# pollen.whit.focal <- matrix(nrow=length(pollen.ages), ncol= ncol(pollen.whit[[i]]$counts))
# for (i in 1:length(pollen.ages)){
#   sampleRow <- which.min(abs(focalTime - pollen.ages[[i]]$age))
#   pollen.whit.focal[i,] <- pollen.whit[[i]]$counts[sampleRow,]
# }
#
# #Stratiplot(pollen.whit[[1]]$counts, pollen.whit[[1]]$chronologies[[1]]$age, ylab="Age (cal BP)", sort="wa", c("poly","g"))
#
# whit.counts <- lapply(pollen.whit, '[[', 4)
# whit.S.obs <- lapply(whit.counts, specnumber)
# whit.shannon <- lapply(whit.counts, diversity, "shannon")
ggiNEXT(iNextRichness, type=1, facet.var="order", color.var="site")
ggiNEXT(iNextRichness, type=1, facet.var="site")
ggiNEXT(iNextRichness)
ggiNEXT(iNextRichness, type=2)
library('vegan')
library('vegan')
library('neotoma')
library('dplyr')
library('sads')
library('iNEXT')
?neotoma
??neotoma
samwell_site <- get_site(sitename="Samwell Cave") #there's one site in the database that matches this
samwell_site
str(samwell_site)
library('vegan')
library('neotoma')
library('dplyr')
library('sads')
library('iNEXT')
samwell_site <- get_site(sitename="Samwell Cave") #there's one site in the database that matches this
samwell_datasets <- get_dataset(samwell_site[[1]]) # this search lists all datasets associated with Samwell Cave. There are 7 different datasets for the site: 4 containing vertebrate faunal information and 3 containing geochronologic information.
samwell_datasets
str(samwell_datasets)
browse.dataset(scpd_dataset)
browse(scpd_dataset)
browse(samwell_datasets)
scpd_dataset <- get_dataset(samwell_site[[1]])$'14262'
scpd_dataset
str(scpd_dataset)
scpd_download <- get_download(scpd_dataset)
str(scpd_download)
scpdFauna <- counts(scpd_download)[[1]]
scpdAges <- ages(scpd_download)[[1]]
head(scpdFauna)
scpdAges
str(scpdFauna)
all(rownames(scpdFauna) == scpdAges$sample.id
)
if (all(rownames(scpdFauna) == scpdAges$sample.id)){
rownames(scpdFauna) <- ages
}
scpdFauna <- counts(scpd_download)[[1]]
scpdAges <- ages(scpd_download)[[1]]
ages <- scpdAges$age
if (all(rownames(scpdFauna) == scpdAges$sample.id)){
rownames(scpdFauna) <- ages
}
head(scpdFauna)
sm_genus <- c("Scapanus", "Microtus", "Neotoma", "Peromyscus", "Sorex", "Spermophilus", "Tamias", "Tamiasciurus", "Thomomys", "Arborimus", "Reithrodontomys", "Sciurus", "Glaucomys", "Myodes", "Aplodontia")
scpdSM <- matrix(nrow=nrow(scpdFauna), ncol=length(sm_genus))
head(scpdSM)
i <- 1
sm_genus[i]
grep(sm_genus[i], colnames(scpdFauna))
scpdFauna[,2]
head(scpdFauna)
for (i in 1:length(sm_genus)){
cols <- grep(sm_genus[i], colnames(scpdFauna))
if (length(cols)==1){
scpdSM[,i] <- scpdFauna[,cols]
}
if (length(cols) > 1){
scpdSM[,i] <- rowSums(scpdFauna[,cols])
}
}
scpdSM
# clean up dataframe
colnames(scpdSM) <- sm_genus
rownames(scpdSM) <- rownames(scpdFauna)
scpdSM
# we can practice building our own functions to count up the species
countFcn <- function(df){
count <- length(which(df>0))
return(count)
}
raw.S <- apply(scpdSM, 1, countFcn)
?apply
raw.S
specnumber(scpdSM)
range(rowSums(scpdSM))
minN <- min(rowSums(scpdSM))
minN
richness <- rarefy(scpdSM, sample=minN, se=T) # because I specified *se=T*, this function outputs both S and se
richness
S <- richness[1,]
S.se <- richness[2,]
plot(S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, floor(max(S+S.se))+1), xlab="Age (cal BP)", ylab=expression('Richness' %+-% 'se'))
# plot the result through time using base R plotting function
plot(S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, floor(max(S+S.se))+1), xlab="Age (cal BP)", ylab=expression('Richness' %+-% 'se'))
#arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
suppressWarnings(arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90))
plot(S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, floor(max(S+S.se))+1), xlab="Age (cal BP)", ylab=expression('Richness' %+-% 'se'))
arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
plot(S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, floor(max(S+S.se))+1), xlab="Age (cal BP)", ylab=expression('Richness' %+-% 'se'))
#arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
suppressWarnings(arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90))
plot(raw.S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, max(raw.S)+1), lty=2, col="gray", xlab="Age (cal BP)", ylab='Richness')
lines(S ~ ages, type="l")
arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
# set colors, from cool for the older (LGM) times and red for the warmer (Holocene) times
ageColors <- rainbow(length(ages), start=0, end=4/6)
rarecurve(scpdSM, col=ageColors)
?estimateR # read more about this function, and how ACE and Chao1 are calculated
e
alt_richness <- estimateR(scpdSM)
alt_richness
# plot them against one another
y.max <- floor(max(richness_all))+1 # set y max for plotting
richness_all <- rbind(alt_richness, richness)
# plot them against one another
y.max <- floor(max(richness_all))+1 # set y max for plotting
mainEstimates <- c('S.obs', 'S.chao1', 'S.ACE', 'S')
plotRows <- match(mainEstimates, rownames(richness_all))
plotCols <- c("black", "blue", "green", "red")
i=1
plot(richness_all[plotRows[i],] ~ ages, type="l", ylim=c(0, y.max), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[i])
for (i in 2:length(plotRows)){
lines(richness_all[plotRows[i],] ~ ages, col=plotCols[i])
}
legend("bottomleft", legend=mainEstimates, lty=1, col=plotCols, bty="n", cex=0.75)
?iNEXT
vignette("Introduction", package="iNEXT")
vignette("Introduction", package="iNEXT")
?iNEXT
iNextRichness <- iNEXT(t(scpdSM))
iNextRichness
D_minN <- estimateD(t(scpdSM), datatype="abundance", base="size", level=NULL)
D_minN
D_minN[which(D_minN$order == 0),]$qD
richness[1,]
ggiNEXT(iNextRichness)
vignette("diversity-vegan")
?diversity
H <- diversity(scpdSM, index="shannon")
Simp <- diversity(scpdSM, index="simpson")
H
Simp
Shan_even <- H/log(specnumber(scpdSM))
Simp_even <- Simp/log(specnumber(scpdSM))
plotCols <- c("black", "blue",  "green", "red")
plot(H~ages, type="l", ylim=c(0, 2), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
lines(Simp~ages, col=plotCols[2])
lines(Shan_even~ages, col=plotCols[3], lty=2)
lines(Simp_even~ages, col=plotCols[4], lty=2)
legend("bottomleft", legend=c("Shannon", "Simpson", "Evenness_Shannon", "Evenness_Simpson"), lty=c(1,1,2,2), col=plotCols, bty="n", cex=0.75)
vignette('sads_intro')
scpdSM[1,]
scpdSM[14,]
scpd1.oc <- octav(scpdSM[1,])
scpd14.oc <- octav(scpdSM[14,])
scpd1.oc
scpd14.oc
scpd1.oc
scpd14.oc
scpd1.oc
scpd14.oc
par(mfrow=c(4,4), mar=c(3,3,2,0)+0.1)
for (i in 1:nrow(scpdSM)){
plot(rad(scpdSM[i,]), prop=F, xlab="", ylab="", xlim=c(1, 13), type="b", pch=16)
mtext(side=2, line=2, "Number of Indiv", cex=0.5)
mtext(side=3, line=0.5, paste0("Time=", ages[i]), cex=0.75)
mtext(side=1, line=2, "Rank", cex=0.5)
}
?rad
scpd1.geom <- fitsad(scpdSM[1,scpdSM[1,]>0], 'geom')
scpd1.ls <- fitsad(scpdSM[1,scpdSM[1,]>0], 'ls')
scpd1.lnorm <- fitsad(scpdSM[1,scpdSM[1,]>0], 'lnorm')
AICtab(scpd1.geom, scpd1.ls, scpd1.lnorm, base=TRUE)
